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SUMMARY 


A numerical method has been used to calculate steady, three-dimen- 
sional, viscous compressible flow fields about slender bodies at angle 
of attack. The method is based on the equivalence principle, and on a 
procedure developed earlier for solving the equations of two-dimensional 
transient continuum motion (applied here to the Navier-Stokes equations) . 
The numerical technique, and the basic features of the equivalence 
principle as applied to case of present concern, are both reviewed herein, 
and the results of calculations of two slender body flow fields are pre- 
sented. 

The equivalence principle relates the axial coordinate of a steady 
three-dimensional flow field to the time coordinate of two-dimensional 
transient flow; thus, a given three-dimensional problem becomes directly 
analogous to a two-dimensional problem of flow about a cylinder (not 
necessarily circular) that deforms with time. It is basic to the 
validity of the equivalence principle that axial velocities do not differ 
significantly at any point from the axial component of the free stream 
velocity, a requirement that precludes the use of a no-slip condition at 
the body surface in the axial direction. Accordingly, the flow fields 
considered in the program, including separation and body vortex develop- 
ment, were computed with a model that permits viscous crossflow together 
with inviscid axial flow. An analysis of the errors introduced by such 
a treatment is presented herein. 

Numerical calculations were made and compared with experimental 

results for an ogive-cylinder and an airplane fuselage configuration. 

6 o 

Flow conditions were = 1.98, = 4.68 x 10 / ft, and a = 10 for the 

6 o 

ogive-cylinder, and M ra = 2.50, R ro = 9.1 x 10 /ft, and a = 15 for the 
fuselage configuration. Results are presented as static pressure distri- 
butions on the body surface, as velocity vector plots, and as contour maps 
of the flow field at selected axial locations. Good agreement between 
theory and experiment was obtained; maximum deviations of numerical 


li 


I 


surface pressures from corresponding experimental values were no more 
than 6 percent of free stream dynamic pressure for both problems. 
However, some differences were noted. Boundary- layer separation and 
body vortex positions differed from experimental locations on the 
ogive-cylinder, and the shock induced by the fuselage canopy was 
predicted at a slightly different location. These differences are 
considered attributable to neglect of axial viscous effects, exclusion 
of turbulence phenomena, and the approximations introduced by the 
equivalence principle in describing inviscid axial flow. 
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1.0 INTRODUCTION 


For highly maneuverable advanced aircraft (cf . figure 1) , the 
prediction of flow-field characteristics for nonsimple geometries at 
high angles of attack and high flight Mach numbers transcends the 
capabilities of linear theories. For example, there is a need for an 
analytical or numerical method that provides an accurate , detailed descrip- 
tion of the vehicle flow field in order to determine both fuselage con- 
figurations and inlet-airframe integration effects. 

Prediction of the three-dimensional flow field about a slender body 
immersed in a supersonic airstream presents a formidable mathematical 
problem to which no exact solution has been obtained, either analytically 
or numerically. Successful numerical methods using the method of charac- 
teristics and the equivalence principle have been developed for three- 
dimensional, steady, inviscid flow fields. The work of Gallo and Rakich 
(ref. 1) is representative of the characteristics approach, while 
Van Dyke's extension '"ref. 2) of the theory introduced by Hayes (ref. 3) 
is representative of an approach that uses the equivalence principle. 

In this research program the equivalence principle was used to 
convert time-dependent viscous, compressible flow calculations in two 
space dimensions to steady motion of a viscous, compressible fluid in 
three space dimensions. The two-dimensional flow calculations were 
obtained by numerically solving the time-dependent compressible Navier- 
Stokes equations. 
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2.0 SYMBOLS 


a 

a 


A 



B 2 (z') 


B'(z') 

B 5 

B(x,y,z) 


C 


C 


D 


C 


L 


C 


P 


c 

P 


c 

V 


d 


d 


Lagrangian coordinate 
Radius 

Surface area 

Abscissa of center of smaller circle of fuselage cross section 
Abscissa of center of canopy circle 

Ordinate of center of smaller circle of fuselage cross section 

Ordinate of axis of the body with respect to the horizontal 
reference line 

Slope of axis of the body with respect to horizontal reference line 
Ordinate of center of canopy circle 
Surface of the body 


Local sound speed 

Drag coefficient, dra g force 

c 


Lift coefficient, 


Pressure coefficient. 


lift force 
q A 

CO C 

P - P 


Specific heat at constant pressure 
Specific heat at constant volume 
Diameter 

Displacement vector 
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E 



m 

M 


M 

| AM ] 


N 

N 

P 

P 

s 

|AP| 


q 

R 

R, 


R„ 


R( j 5 k) 
R'(j,k) 


Specific internal energy 

Bounded oscillating function used to determine discretization error 
Stagnation enthalpy 

Direction cosines of the normal to the body's surface 
in the x,y, and z directions, respectively 

Mass 

Mach number 

Momentum vector 

Absolute error in Mach number 

Direction cosines of the normal to a shock surface in the x, y 
and z directions, respectively 

Cycle number 

Mesh point density 

Pressure 

Stagnation pressure 
Absolute error in pressure 
Pitot pressure 

2 

Dynamic pressure, q = (1/2) pU 
Reynolds number 

Radius of smaller circle of fuselage cross section 
Radius of larger arc of fuselage cross section 
Radius of canopy circle 

Position vector of the mesh point about a fuselage cross section 
Position vector of a mesh point about a circular cross section 
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s 

S(x,y,t) 

t 

At 

T 

u 

U 

U 

v 

V 

V 
w 
W 

x 

x* 

X 

y 

y* 

z 

z' 

a 


Entropy 

Surface describing system of shock waves about a body 
Time 

Time step 
Temperature 

Velocity component in x direction 
One-dimensional velocity, U(a,t) 

Velocity vector 

Velocity component in the y direction 

Specific volume 

Volume 

Velocity component in the z direction 
Work rate 

Coordinate normal to plane consisting of the 7%° reference line 
and the perpendicular to that line 

Coordinate normal to the plane consisting of the horizontal 
reference line and the perpendicular to that line 

Lagrangian coordinate, X = X(a,t) 

Coordinate normal to 7%° reference line 

Coordinate normal to the horizontal reference line 

Coordinate along the 7%° reference line 

Coordinate along the horizontal reference line 

Free-stream angle of attack with respect to body axis 
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Local angle of attack with respect to a horizontal reference line. 



/ velocity in y* direction 
^velocity in z‘ direction 


Sideslip angle, 


tan 


( velocity in x* direction 
velocity in z 1 direction 


c 

Ratio of specific heats, — ^ 

c v 


Discretization error in a property (eq. (29)) 


Density 

Maximum slope of surface of body with respect to the free-stream 
flow direction 


Stagnation condition 

Property downstream of a shock 

Free-stream condition 

Lagrangian coordinate at time zero 

Reduced variable 

Actual variable 


Variable associated with body cross section 



3.0 THE NUMERICAL METHOD 


3.1 FORMULATION OF EQUIVALENCE PRINCIPLE 

The two-dimensional time-dependent equations of motion (independent 
variables x,y,t) can be transformed directly to steady-state three- 
dimensional (x,y,z) equations for slender bodies through the Hayes equi- 
valence principle. The equivalence principle relates the steady-state 
flow field over a slender body to an equivalent time-dependent flow field 
in one less space dimension. In particular, steady three-dimensional flow 
about a body can be reduced under certain conditions to two-dimensional 
time-dependent flow in a plane normal to the free-stream flow direction. 

A simple example of the application of the equivalence principle is 

illustrated in figure 2, namely that of steady flow over an axisymmetric 

body at zero angle of attack. The steady-state two-dimensional axially 

symmetric flow field is analogous to time-dependent one-dimensional 

axially symmetric flow in a plane caused by an expanding cylinder. The 

cross-section plane of unsteady analogy moves downstream with free-stream 

velocity and the outline of the moving boundary is given by the trace of 

the original shape in the cross-section plane. At the station z (see 

fig. 2) , the steady-state flow field in the cross-section plane is analogous 

to the flow field about the expanding cylinder at the time, t = z/U m . To 

make the equivalence analogy valid, the normal velocity of the expanding 

cross section U , is given by the product of the free-stream velocity and 
n 

the tangent to the surface. Far from the body, free-stream conditions 
are imposed (i.e., P = P ra , p = p ro , u = v = 0). 

According to convention, the free-stream flow direction is chosen 
as the principal axis along which the linearization approximations for 
equivalence are made. However, for bodies at angle of attack as illus- 
trated in figure 3, it is more convenient to choose the body axis as the 
principal axis. In planes normal to the body axis, the cross sections 
can usually be defined in terms of a few simple analytic expressions, 
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which result in simple and accurate surface boundary conditions, and 
facilitate the development of finite difference meshes about these cross 
sections. Because of this choice, free-stream crossflow will exist as an 
upstream boundary condition, but the choice of axis in no way affects the 
accuracy of the results. This point is further discussed in Section 3.3. 

It is convenient to discuss the equivalence principle in terms of 
flow about a slender body of revolution at angle of attack, O'. Let the z 
axis coincide with the axis of the body, and let the cross section of the 
body lie in the x,y plane (see fig. 3). It is also assumed that the z 
component of the local velocity vector is approximately equal to the z 
component of the free-stream velocity vector, that is, 

w = U cos a (1) 

00 

Under this hypothesis, it follows that time-dependent flow in the x,y 
plane transforms to time-independent motion in x,y,z space according to 
the equations 


z = cos at (2) 

y = y (3) 

x = x (4) 

This result — the equivalence principle — is derived for an inviscid fluid 
in Appendix A. 

In order for equations (2), (3), and (4) to represent a valid map- 
ping between a flow field in x,y,t space and a flow field in x,y,z space, 
the time-dependent solution in the x,y plane must be augmented by a three- 
dimensional boundary condition at the body surface 

1 U cos a+Xv+£u = 0 (5) 

z ta y x 
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where i 


1 , i are the direction cosines of the normal to the body 
y z 


surface in the x,y, and z directions, respectively. Equation (5) states 

that at the body, the component of the local velocity normal to the body 

surface is zero. The application of the boundary condition ^ equation (5)j, 

in the x,y plane implies that the body cross section varies in time and 

that the velocity normal to the surface is equal to -i U cos a/ / j£ 2 + j£ 2 

z oo v y x 

a normal velocity expression that reduces to the product of the axial 
velocity (U M cos a) and the local tangent to the surface. 


A second boundary condition, a constraint on tangential flow at the 
body surface, must also be specified. Since the full time-dependent Navier- 
Stokes equations are solved in the x,y plane, a no-slip boundary condition 
is imposed at the cross section of the body. We require that 


vi - ui =0 (6) 

x y ' 

It is seen in Appendix A that the boundary equations ^equations (5) and (6)j 
are compatible with the assumptions inherent in the equivalence principle 
and permit a time-dependent viscous calculation to be made in the x,y plane. 
Thus, with the inclusion of boundary equations (5) and (6), a solution of 
the two-dimensional time-dependent Navier-Stokes equations is also the 
solution to a problem of steady three-dimensional viscous flow. 

The three-dimensional boundary layer on the body surface has been 
replaced by a two-dimensional boundary layer since the no-slip condition is 
imposed at the body surface only in the crossflow direction. This 
ad hoc description of the flow field is poor in the boundary layer. 

However, its justification lies in the fact that it provides a mechanism 
for effecting flow separation and the subsequent development of spiral 
vortex sheets on the lee side of a body at angle of attack. To be sure, 
solving the compressible Navier-Stokes equations in a crossflow plane is 
much more realistic than the conventional approach of solving the Laplace 
equation in this plane with experimentally (or otherwise) determined 
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circulation strengths. However, the accuracy of the viscous crossflow 
description must he carefully documented by comparing the numerical 
results with other calculations and with experimental data. 

The computational procedure is then as follows. Initially we have 
a uniform flow field in the x,y plane with a velocity in the y direction 
equal in magnitude to sin The slender-body cross section grows with 
time in the x,y plane. At discrete times solutions of the full time- 
dependent Navier-Stokes equations are obtained such that the cross-section 
surface boundary conditions ^equations (5) and (6)j are satisfied; also, 
the flow far from the cross section is uniform with a crossflow velocity 
in the y direction equal to sin a. The time-dependent solution in the 
x,y plane is then related to the steady three-dimensional flow field 
about the body through the transformation equations (2), (3), and (4). 

3.2 DESCRIPTION OF THE NUMERICAL METHOD 

The numerical method used to solve the time -dependent Navier-Stokes 
equations in two space dimensions is embodied in a computer code called 
"AFTON 2PE' r . The finite difference equations in the AFTON 2PE computer 
code are based on a physical model of the continuum. The need for such a 
model as a basis for numerical calculations of continuum motion has long 
been recognized (refs. 4,5 and 6) } and was filled after considerable effort 
for continuum motion in all its generality, the model for two-dimensional 
flow actually preceded the general case, and is recalled in some detail 
below. The two-dimensional continuum model is not specific to the par- 
ticular phenomenon under study here, nor even to the entire class of 
motions governed by the Navier-Stokes equations; once a physical model of 
the continuum is verified, its range of application extends to many 
flow situations. 

In the conventional use of physical models, it is assumed that the 
partial differential form of the Navier-Stokes equations applies in general. 
A physical model is then postulated for the particular phenomenon to be 
studied and the Navier-Stokes equations are reduced to ordinary differen- 
tial equations from which we obtain either an analytical or numerical 
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solution. The Lees-Reeves (ref. 7) near-wake model is a good example of 
a conventional physical model. This model is concerned with the base 
flow and near-wake subregions shown in figure 4. The model recognizes a 
constant pressure shear-layer near the body separation point and a vis- 
cous inviscid flow interaction near the wake-stagnation point. The bound- 
ary layer equations are assumed to govern the flow in both regions. The 
flow fields in both regions are solved separately and then joined together 
through the pressure gradient in the streamwise direction. 

For illustrative purposes, the space-time continuum model employed 
in this research is presented in detail for the case of time-dependent one- 
dimensional flow of a compressible, inviscid fluid. This model was implied 
by a set of finite difference equations developed by von Neumann and Richt- 
myer (refs. 8 and 9) . A careful analysis of the von Neumann-Richtmyer 
equations led to the continuum model employed here (refs. 4,5 and 10). 

Finite difference analogs of the continuity equation, momentum 
equation, and the First Law of Thermodynamics will be derived in a Lagrangian 
coordinate system. Let a be the Lagrangian coordinate, and X(a,t) be the 
Euler ian coordinate. That is, X(a,t) gives the position at time t, of a 
fluid element that was originally at position a. Consider the Lagrangian 
coordinates, a x , a 2 , and a^ [a^ = % (a x + ag)] shown in figure 5. 

Since the system is Lagrangian, the mass between the trajectories labeled 
a-!^ and a 3 remains a constant. Let the one-dimensional space continuum be 
represented by a discrete set of zones, designated "thermodynamic" zones. 

At time zero, let the boundaries of these zones be spaced at a constant 
interval Aa along the X axis [see fig. 5 3 and denoted by a^(^=0, 1,2,3, 

...,L). Let each of the thermodynamic zones be one unit high and one 
unit wide. Consider another set of zones designated "momentum" zones 
superimposed on the thermodynamic zones, such that each momentum zone 
surface always divides the mass of the thermodynamic zone which contains 
this surface in half. A schematic diagram showing a momentum zone and 
two thermodynamic zones at the time t is shown in figure 5. The 


10 



Eulerian coordinates which define the two thermodynamic zones are 
X(a^ ^,t), X(a^, t), and X(a^_^,t), while the Eulerian coordinates of the 
momentum zone are X(a ,,t) and X(a . , t) . 

ju"**2 X* i 'Z 


The thermodynamic variables, such as internal energy, specific vol- 
ume, and pressure, are assumed constant throughout a thermodynamic zone. 

The velocity is assumed constant throughout a momentum zone. Thus, thermo- 
dynamic variables are in effect centered at the momentum zone surfaces; 
that is, P^ = P(a^ i„,t) denotes the pressure in the thermodynamic zone 
X(a 1 ,t), X(a ,t) which effectively acts at X(a , ,t) [see fig. 5 ]. 
Momentum zone variables are in effect centered at the thermodynamic zone 
surfaces; that is, U = U(a ,t) denotes the particle velocity in the mo- 

Xj Xj 

mentum zone X(a^ j.,t), X(a^ ^,t) which effectively acts at X(a^, t). 


Since an explicit formulation is the goal, the variables of motion 
are not only displaced spatially (as discussed above), but they are dis- 
placed in time as well. Let At denote the uniform time interval and 
t n (n=0,l,2, . . . ,N) denote the time after n uniform time intervals. The 
variables associated with thermodynamic zones are defined at integer times; 
that is, = P(a^ denotes the pressure at the time t n . The mo- 

mentum zone variables are defined at half-integer times; that is, 


U 


n-% 




n-% 


= U(a t 2 ) denotes the particle velocity at a time t 2 , and the 


Eulerian coordinate position X(a , t z ) . 

Xj 


Finite difference analogs of the continuity, momentum, and first law 
equations follow directly from the physical model of the continuum pre- 
sented in figure 5. Since a Lagrangian coordinate system is employed, 

the zones of figure 5 will be displaced continuously from their initial 
positions; that is, a 1 ,...,a^ to coordinate positions X(a x , t) , . . . ,X(a^ , t) . 
Let us calculate the properties at the time, t n , for the thermodynamic 
zone having Lagrangian coordinates a^_-^,a^ and the momentum zone having 
Lagrangian coordinates a , ,a , . The initial values for this calculation 

^ Xj^'z xsr-% 

are 

n-l _n-l n-1 p n-l n-1 n-1 

je-v v SL-k’ AMs’ z+k’ v n+h 
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for the thermodynamic zones and 

x* n , x* x* u“-J, C % , 

i-V j0+l’ .0-1’ A ’ jfrKL 

for the momentum zones. The objective is to update these variables by one 
time step in an explicit manner. 

The finite difference analog to the conservation of mass is derived 
from the expression for the volume change of thermodynamic zone a , a . 

Xj 1 Xj 


V* a - vT, 1 = (u”"* - U*iV 

j0-% j0-% \ JL -0-1 / 


(7) 


where 


V 


n-1 

- 0 -% 


volume of zone a. , ,a„ at time t 
- 0 - 1 ’ -0 


n-1 


V 


, volume of zone a .,a at time t 

X/““ 1 Xj 

Based on equation (7), and the fact that the mass of material in the 

thermodynamic zone a 1 ,a has the constant value p Aa, we find that 

Xj** IX/ O 


t n n-l N 

v j0-% ~ v -0-% 

At 


ur % 


n-% 

- 0-1 


Aa 


( 8 ) 


where p Q is the density at time zero, and v^ ^ is the specific volume of 
material in the zone at time t n (note that (p Aa)v” , = ,) . 

O X '""'2 X )“^2 

For an inviscid, adiabatic fluid, the first law equation for a sys- 
tem in equilibrium is applicable: 


DE Dv 

Dt V Dt 


(9a) 


The First Law finite difference analog to equation (9a) is derived by 
first writing down a finite difference analog to the term -P(Dv/Dt) and 
then equating this to the rate of change of internal energy in the zone. 
Since the zone mass is constant and the pressure is homogeneous in a 
thermodynamic zone, the term -P(Dv/Dt) becomes 
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where 


= p o Aa 

On the basis of equation (9b) , the finite difference First Law equation 
becomes 


’o( E i 




,n-l 


where is the specific internal energy of thermodynamic zone a^ ^,a^. 

The equation of state for a perfect gas is: 

Pv 

E = (1] 

Finally, the finite difference equation of motion for momentum zome 
a _e-%> becomes: 


M n n-1 
M . - M „ 

A l 


P n “% _ p n-% 
A-% i+k 


P n "% = l/ P n + P n-l\ 

l-h 2\ i-h l-h! 

P n “% = l/ p n . p n-l\ 

£+% 2\i+k l+\) 


and defines the momentum in zone a at time t n . 

Xj Xj 

Equations (8), (10), and (12), derived from the physical model pre- 
sented in figure 5 are solved in the following manner. First, the 

specific volume, v^, is calculated from equation (8) . Then, the 
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are calculated from 


pressure and specific internal energy E^ 

equations (10) and (11). The continuity and first law equations 
are then solved for specific volume and pressure in all thermodynamic 
zones. Based on the pressure field, equation (12) can be solved for the 
momentum . The particle velocity U 1 ^" 1 " 2 is found from the forward 
extrapolation formula 


, , m : a 

U n+% = 2-1 _ u"-% 

1 1 


(13) 


where 


m „ = o Aa 

Integration of the particle velocity U 2 with time establishes the new 
Eulerian position, X . This can be done for all momentum zones. Thus, 
all the initial values cited above can be updated by one time step and 
the process repeated. 

The success of this method coupled with the failure of other more 
recent numerical schemes to improve upon it, led Trulio and Trigger to a 
careful analysis of these finite difference equations in order to dis- 
cover the reasons for their success. They found that the finite differ- 
ence equations possessed the same self-consistency property of form as 
the original differential equations from which they were derived, that is, 
the finite difference equations for momentum conservation and the first 
law implied an exact conservation of total energy (internal and kinetic) 
finite difference equation (ref. 4). In other words, the continuity, 
momentum, first law, and conservation of total energy relations are 
redundant by one. In most other numerical schemes if one tries to derive 
a conservation of total energy relation from finite difference analog of 
momentum conservation and the first law, error terms, usually assumed to 
be "second order", are produced. These numerical error terms are believed 
to be a primary source of the difficulties encountered in many other 


14 



numerical schemes. Appendix B demonstrates the self-consistency property 
of the finite difference equations (8) to (13). 

In addition to self-consistency of form, the numerical method has 
the following properties: 

(a) Since the method is explicit, stability criteria must be met 
in order to obtain physically meaningful numerical results. 

In general, the time step must be small enough so that a 

sound signal cannot cross a zone in a time step (refs. 8 and 9). 

(b) The numerical error has been correlated for this scheme. It 
has been found that the absolute error in a property is in- 
versely proportional to the linear mesh point density to the 
three-halves power (ref. 10) . 

The physical model, from which equations (8) to (13) have been 
derived has been extended to two space dimensions (ref. 11) . In that work, 
specific finite difference equations were formulated and their self- 
consistency properties demonstrated. Because they provide a base for the 
present calculation, a brief description of the derivation of these equa- 
tions will be presented. 

Consider two-dimensional flow of a viscous, compressible fluid in 
the x, y plane. As in the one-dimensional case, let us divide the continuum 
into two types of zones; namely, quadrilateral zones and momentum zones 
(see figure 6). The assumptions governing this analysis are as follows: 

(a) All zones are polygons. 

(b) The density, specific internal energy, pressure, stress tensor, 
and velocity derivatives are homogeneous in a quadrilateral 
zone. 

(c) The velocity vector is homogeneous in a momentum zone. 

(d) The zones have unit thickness normal to the plane of motion. 
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As in the one-dimensional case the continuity and first law equa- 
tions are solved for each quadrilateral zone to determine the density, 
specific internal energy, and stress. When the stress tensor depends 
only on strain, the first law equation corresponds to that of thermo- 
dynamics. However, for a stress-rate of strain dependence, the first law 
equation is a relationship involving internal energy which results from 
subtraction of the kinetic energy equation from the conservation of total 
energy equation. Based on the stress tensors in each of the four quadri- 
lateral zones which comprise a momentum zone (see fig. 6), the momentum 
equation is solved for each momentum zone. To accommodate the equivalence 
principle boundary conditions, the finite difference equations in 
AFTOK 21’ E are written in a generalized coordinate system where the four 
mesh points comprising a quadrilateral zone can move with arbitrary 
velocity. The numerical method is described in detail for a generalized 
coordinate system and an Eulerian coordinate system in two places 
(refs. 11 and 12). 

All of the salient properties of the numerical method described in 
the one-dimensional example are preserved in two-dimensions and will also 
be preserved in three spatial dimensions; in that respect the numerical 
theory is internally consistent. The momentum and first law finite dif- 
ference equations in two dimensions imply an exact finite difference 
equation for total energy. The stability criteria are the same as in the 
one -dimensional case. The one-dimensional correlation of absolute 
numerical error in a property extends directly to two dimensions. In two 
dimensions the absolute numerical error is inversely proportional to the 
area mesh point density to the three-quarters power (ref. 12). The 
AFTON 2PE finite difference equations at interior mesh points are pre- 
sented in Appendix C. 

3.3 ORDER OF ERROR OF THE NUMERICAL METHOD 

In this section the order of error of the numerical method is 
determined. The numerical error is defined as the absolute difference 
between the actual value of a quantity and the value computed numerically 


The specific internal energy is the difference between total and 
kinetic energy. 
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The error in this method stems from three sources, the equivalence prin- 
ciple, discretization errors, and neglect of axial viscous effects and 
turbulence phenomena. The errors inherent in the equivalence principle 
lie in the neglect of the velocity perturbation along the free-stream 
flow direction and have been deduced by Van Dyke (ref. 2). The discreti- 
zation error is defined as the difference between the exact solution to 
the continuum motion equivalence principle equations for a given system, 
and the values of the flow variables computed from a finite difference 
approximation to these equations of motion; discretization error results 
basically from the substitution of a discrete set of points for the 
space-time continuum. The errors introduced by neglecting axial viscous 
effects and turbulence phenomena manifest themselves in the accuracy with 
which separation point locations and vortex center trajectories can be 
predicted. Since the determination of the accuracy of the method is of 
primary importance, all three sources of error will be discussed at some 
leng th. 

3.3.1 Equivalence Principle Errors 

To deduce the order of error in the neglect of the velocity pertur- 
bation along the free stream flow direction. Van Dyke introduced what he 
termed "reduced" independent and dependent variables, and he redefined 
the functions describing the body and shock-wave surfaces. To obtain 
these reduced variables, consider a coordinate syscem where the z axis is 
aligned with free-stream flow direction. Let the surface of the body be 
described by B(x,y,z) = 0, and let the complete system of shock waves be 
described by S(x,y,z) = 0. In this coordinate system the reduced 
variables are as follows: 



( 14 ) 
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u = 

V = 

w = 
p = 

p = 

B = 

S = 

where t is the maximum slope of the body surface with respect to the free- 
stream flow direction. This transformation of variables was introduced 
into the continuity, momentum, and first law equations for an inviscid 
fluid. Reduced parameters were considered of order one or less and terms 
that contained j 2 explicitly were discarded. A set of reduced equations, 
linear in w and derivatives with respect to z", resulted. 

If the reduced equations are rewritten in terms of the actual un- 
barred parameters, and the substitution z = U^t is made, the time-depend- 
ent equations of motion in two space dimensions result. 

The boundary conditions for the reduced equations are as follows: 

At the body surface 

at B = 0 (19) 


( 20 ) 

The parameters M and r of the full problem enter the reduced 
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problem only in the combination M^t, which appears only in the upstream 


= o 

Ox dy dz 


Far upstream of the body 

P - 


YM 2 t 3 


as 


U m ru(x,y,z) \ 


U Q Tv(x,y,z) > 

(15) 

U m [l + T 2 w(x,y,z)] ) 


P vM 2 T 2 P(x,y,z) 

co' oo 

(16) 

pj>(£,y,2) 

(17) 

B(x,y,z) ) 



(18) 

S(x,y,z) ) 



18 


boundary condition on P (eq. 20) . Since P must be of order one or less 
for these equations to be valid, M^T, must be of order one or greater for 
this theory to be consistent far upstream of the body. If this consistency 
condition is satisfied, the maximum error in a reduced variable must be 
of the order t 3 , because terms of order t 3 have been omitted in the 
reduction. Since M t is of order unity, an error of order t 3 implies an 
error of order 1/(M 3 ). Hence, the theory becomes more accurate as the 
free stream Mach number goes up. Since supersonic flow problems are 
solved in this paper, this source of error, although small, will not be 
negligible. Based on the t 3 error law, the maximum errors in the static 
pressures and local Mach numbers are derived below. 

The absolute difference between the actual reduced pressure, P , 
and the reduced pressure calculated from this theory, P, is 

|P a - p| w t 3 (21) 

Combining equations (16) and (21) yields 




T 


3 


where q^ is the free-stream dynamic pressure [(%)p Q U®]. Therefore, the 
absolute error in the pressure is as follows: 



( 22 ) 


Since experimental and numerical local Mach numbers are compared in 
this research effort, it is important to determine the order of error in 
the Mach number. The Mach number, M, is defined as 


M 


(u 3 + v 3 + w s ) % 
C 


(23) 
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where C is the local sound speed. From Bernoulli's equation, the local 
sound speed can be expressed in terms of the components of local velocity 

C 2 = C 2 - Y -" - " |u 3 + v 3 + w 3 - U 2 ] (24) 

where is the free-stream sound speed. Introducing the reduced veloci- 
ties u, v, and w from equations (15) into equations (23) and (24) yields 
the relation for the Mach number in terms of reduced velocities. 


2 


M = 


M 


-2 


(1 + t 4 w + 2-r a w + t 3 v 2 + t 3 u 3 ) 

' ( X T i ') + 2w + v 2 + u 2 )J 2 


(25) 


Since u, v, and w are of order one, terms in t 4 could be neglected in 
equation (25) . 
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(26) 


Equation (26) can be linearized in terms of changes in u, v, and w 
by employing a first-order Taylor’s expansion. Let M represent the 

3 . 

actual local Mach number and M represent the value computed by this numeri- 
cal method. From a first-order Taylor's expansion 
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The coefficients of equation (27) can be determined by differentiation of 
equation (26) and evaluation of the derivatives at the values u = v — 1, 
w = 0, since barred planar quantities are assumed of order one, and the 
perturbation velocity w is assumed to be near zero. The final equation is as 
follows: 
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where 


M = |M a - M| 

3.3.2 Discretization Errors 

The discretization error in the two-dimensional time-dependent 
equations of Section 3.3.1 is based on a general rule that relates the 
error to the density of mesh points employed in the numerical integration 
process. The derivation of this error rule can be found in references 
10 and 12; the relation is as follows: 

-3/4 

e = f 2 « (29) 


where e is the absolute discretization error in a given property, N is the 
two-dimensional mesh point density (i.e., number of zones per unit area of 
the x,y plane) , and f 2 is a bounded oscillatory function. If the same 
problem is run with two different meshes, that is, a medium and fine mesh; 
the function f 3 can be evaluated from the values of a property and mesh 
point densities of the medium and fine meshes at the same point in space 
and time. For the case where the static pressure error is required. 


^3 




(30) 


where P. is the fine mesh pressure, P is the medium mesh pressure, N is 
the mesh point density of the medium mesh, and is the mesh point density 
for the fine mesh. 

3.3.3 Errors Due to Neglect of Axial Viscous Effects and Turbulence 

The error introduced by neglecting axial viscous effects and turbu- 
lence phenomena is the most difficult of the three sources of error to 
evaluate. Therefore, no formal attempt was made to evaluate this error 
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on analytical and/or numerical grounds. The numerical results for the 
problems solved in this paper were compared to experimental data. 
Deviations in separation point locations and vortex center trajectories 
were noted, and on the basis of physical arguments, the primary causes of 
the deviations were explained. 

In this research effort the effects of mesh point density on the 
discretization error were not investigated; each problem of the program 
was run with only one mesh. Therefore, the order of error of the numeri- 
cal method, which was established by comparing the numerical results with 
experimental ones, included all three sources of errors. 
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4.0 CALCULATIONS MADE AND RESULTS OBTAINED 

4.1 DESCRIPTION OF PROBLEMS SOLVED 

The AFT0N-2PE computer code, modified to accommodate the equivalence 
principle boundary conditions, was applied first to the flow field about 
an ogive- cylinder configuration and second, to a fuselage geometry repre- 
sentative of an advanced tactical fighter plane. In both problems, air, 
represented as a gamma law gas ( y = 1.4), was considered and adiabatic flow 
was assumed throughout the flow field. For the ogive-cylinder problem 
the free-stream Mach number was 1.98, the angle of attack was 10° with 
respect to the axis of the body, and the free-stream Reynolds number was 
4.68 x 10 /ft. For the fuselage problem the free-stream Mach number was 
2.5, the angle of attack was 15° with respect to the horizontal, and the 
free-stream Reynolds number was 9.1 x 10 /ft. In this section the cross- 
sectional shapes for both problems are defined and in the next section the 
finite difference meshes generated about these cross sections are described. 

The axisymmetric ogive -cylinder configuration is composed of an 
ogive which is three maximum cylinder diameters long and a cylinder 7.3 
diameters long. This configuration is schematically illustrated in 
figure 7, where the equations which describe the variation of the radius 
of the body with axial distance are also indicated. In this problem the 
axis of the body, which is straight in this case, was chosen as the 
principal axis for the equivalence analogy between the steady and unsteady 
flows . 

The fuselage had a drooped nose which resulted in a curved 
central axis of the body. This geometry is schematically illustrated 
in figure 8. The central axis of the body is composed of a straight 
portion, inclined 7%° from the horizontal, and a curved portion which 
begins at a horizontal station of 15 inches (see figure 8.) To avoid 
introducing curvature effects into the equations of motion, the 7%° 
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reference line was chosen as the principal axis for the equivalence 
analogy. In this coordinate system, the angle of attack of 15° with 
respect to the horizontal becomes 7%° with respect to the 7%° reference 
line; hence, 02 = 7%°. The canopy, included in this problem, is 
indicated in the cross sections normal to the 7%° reference line shown 
in figure 8. 

The fuselage configuration has cross-sectional shapes normal 
to the 7%° reference axis whose peripheries can be approximated by 
circular arcs and straight lines. In fact, the fuselage cross section 
is circular to a horizontal station of 9.44 inches from its nose, is 
asymmetric between horizontal stations 9.44 and 10.88, and includes a 
circular canopy between horizontal stations 10.88 and 30.68. The fuse- 
lage cross section at a horizontal station at the canopy location is 
shown schematically in figure 9. The parameters describing this peri- 
phery are also indicated in the figure. These parameters have been 
curve-fitted as functions of distance along the central axis of the 
fuselage. The curve fits of the cross-sectional parameters are pre- 
sented in Appendix D. 

4.2 MESHES USED 

A subroutine of the AFTON 2PE computer code has been developed 
for generating finite-difference meshes around an asymmetric half- 
body of a general fuselage-shaped cross section. This subroutine is 
based on previous work on finite-difference mesh development for a 
circular cylinder (ref. 12) . The general cross-sectional shape of the 
half-body is assumed to consist of two circular arcs, two straight- 
line segments tangent to the circular arcs, and a circular canopy 
(see figure 9) . The procedure adopted in the calculation was as follows: 
The area of the half-body was computed and a half-circle of equivalent 
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area was located with its center at the coordinate origin. The finite- 
difference mesh for this half-circle was calculated from a modified 
stream function and potential function from potential flow theory 
about a cylinder. The problem was to transform the cylinder mesh into 
a new mesh around the asymmetric half-body according to some suitable 
rule of transformation. In the subroutine developed, each half-circle 
mesh point (designated hereafter as an "unprimed" mesh point) was 
transformed into a half-body mesh point (designated hereafter as a 
"primed" mesh point) in the following manner. First, the periphery 
of the half-body shape was divided into as many equal arcs as the 
half-circle periphery. Then, surface vector displacements were 
obtained between corresponding unprimed points on the half-circle 
and primed points on the asymmetric half-body. Based on these surface 
vector displacements, unprimed mesh points in the flow field were 
displaced to their primed locations. Consider an unprimed mesh point 
in the flow field having a position vector R(j,k), where the integer, 
k, is associated with a potential-like line, and the integer, j, is 
associated with a streamline-like line. Let the surface vector 
displacement corresponding to the same k line be :lenoted as d(k) . 

The position vector of the primed point R'(j,k) is determined from 
the following equation: 


R'(j,k) = R(j,k) + i(k)[i^-~] (31) 

where a is the radius of the half-circle and R(j,k) is the magnitude 
of the position vector R(j,k). 


25 


inn 1 1 


The finite-difference meshes for the ogive-cylinder were composed 
of 35 j lines and 90 k lines and continuously deformed with the radius of 
the body. The initial radius was 0.00046875 foot and the maximum radius 
was 0.046875 foot. Since the ogive -cylinder cross section is circular, 
only the equations of the mesh generating subroutine which pertained to 
the circular cylinder mesh points (unprimed mesh points) were used to 
continually calculate new meshes as the radius of the cross section 
changed. The finite-difference mesh in a cross-sectional plane normal to 
the body axis at a radius of 0.01 foot is shown in figure 10. The finite- 
difference mesh corresponding to the maximum radius (0.046875 ft) is 
indicated in figure 11. 

The finite-difference meshes for the fuselage configuration were 
composed of 35 j lines and 99 k lines and continuously deformed as the 
body cross-sectional shape deformed. The mesh- generating subroutine of 
the AFTON 2PE computer code continually calculated new meshes as the 
fuselage shape deformed. The finite-difference mesh in the cross-sectional 
plane normal to the central axis at a horizontal station 7 inches from the 
nose of the fuselage is shown in figure 12. At this station the fuselage 
cross section is circular. The finite-difference mesh corresponding to a 
horizontal station 25 inches from the fuselage nose is shown in figure 13. 
The canopy is also indicated in this figure. 

4.3 BOUNDARY CONDITIONS AND INITIAL CONDITIONS 

In the x,y planes, the finite-difference meshes are bounded by an 
upstream boundary, a lateral boundary, and a boundary composed of the 
symmetry line of the cross section and the body cross section itself. 

The density and specific internal energy are given their free-stream 
values at the upstream boundary while the velocity of material normal to 
this boundary, v , is evaluated from 

v = U sin a (32) 
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where a is the angle of attack with respect to the principal axis. 

Equations (5) and (6) are satisfied at the body's surface, while the fluid 
is allowed to slide without friction at the system's lateral boundary and 
symmetry lines. The downstream boundary condition was based on the method 
of characteristics. It was used previously in low-speed wake-flow calcula- 
tions and gave a good approximation to the flow in this region. References 
12 and 13 describe this downstream boundary condition in some detail. 

The initial conditions for both problems consisted of a uniform flow 
field in the (x,y) plane with zero x-component the velocity and the y- 
component the velocity given by equation (32) . The density and specific 
internal energy were given their free stream values. In both problems 
the initial cross sectional radius was approximately one-hundredth of the 
maximum cross sectional radius. 

4.4 RESULTS OF OGIVE-CYLINDER PROBLEM 

The ogive-cylinder problem was the first attempted in this study. 

This problem was a good test of the method in that vortices had been 
experimentally observed to occur on the leeward side of the body. The 
crossflow Mach number was 0.344 and the crossflow Reynolds number was 
0.7625 x 10^ based on the maximum diameter of the body. The problem was 
run 3406 cycles (i.e., time steps) on the UNIVAC 1108 computer, requiring 
approximately 10 hours of computer time. Solutions were obtained from 
the body's nose to an axial station 8.35 maximum cylinder diameters down- 
stream. At this cross-sectional plane a well-developed pair of vortices 
were calculated. The problem duplicated the body geometry and free-stream 
conditions of a wind-tunnel test (ref. 14). In general, agreement between 
numerical and experimental data was good, providing evidence that the 
numerical method is applicable to bodies for which separation and subse- 
quent development of spiral vortex sheets occur. 

4.4.1 Qualitative Behavior of Numerical Flow Field 

To investigate the qualitative behavior of the flow field, it was 
convenient to exhibit the data in the form of vector plots of the velocities 
of the fluid particules; the tail of each vector corresponds to a mesh point. 
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The sequence of events as one moves down the axis of the body is described 
by figures 14 to 17. In figure 14 the flow field is shown at a station 
0.502 maximum cylinder diameters from the nose of the body. Although the 
finite-difference mesh is relatively coarse with respect to the radius of 
the body at this station, the bow shock is indicated as well as the expan- 
sion which occurs on the leeward side. In figure 15, at a station 2.99 
maximum body diameters, where the ogive section ends, the bow shock is 
better defined. Figure 16 shows the flow field at a station of 4.92 maxi- 
mum diameters, where separation first appears on the leeward side of the 
body. The spiral vortex sheets that develop on the leeward side of the 
body are indicated in figure 17. The formation of the bow shock, the 
leeward expansion, and subsequent development of the spiral vortex sheets 
are all in qualitative accord with experimental observations. The accur- 
acy of these numerical results is considered below. 

4.4.2 Surface Pressure Comparisons 

Numerical pressure distributions around the body are compared to 
experimental data at various axial stations in figure 18. As can be seen 
from the figure, quantitative agreement is achieved from the body nose to 
an axial station 4.92 diameters. At this station the numerical data indi- 
cate that separation had occurred on the leeward side of the body (see the 
velocity vector plot in fig. 16) , whereas the experimental data showed 
leeward separation at a station approximately 6,00 diameters aft of the 
body's nose. As a result, the numerical pressure data on the leeward 
side of the body differed slightly from the experimental data between 
stations 4.92 and 6.00 diameters down the body axis. At stations greater 
than 6.00 diameters from the cylinder's nose, the experimental data also 
showed separated flow, and the numerical and experimental pressure coeffi- 
cient data were in close agreement 7.63 diameters from the nose of the 
cylinder. 

4.4.3 The Flow Field About the Ogive-Cylinder 

The circumferential positions of the separation points and vortex 
centers vary with axial location in a manner determined experimentally by 
measurements of pitot pressure at various body cross sections. Corres- 
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ponding separation point and vortex center positions were also found from 
velocity plots of the numerical flow field. The calculated separation 
points and vortex centers were found to lie about 20° closer to the wind- 
ward side than the corresponding experimental values as seen in figure 19. 
For example, in a crossflow plane 8.3 diameters down the body axis, the 
numerical separation point and vortex location were 109° and 137°, 
respectively, while the corresponding experimental values were 130° and 
160°. The discrepancy between measured and calculated separation point 
and vortex center positions, while not very large, is perhaps the least 
satisfactory feature of the calculation. The difference is due to 
viscous effects of axial motion being neglected as well as to the exclu- 
sion of turbulence phenomena. 

The flow field predicted numerically is laminar. However, the 
actual flow field is turbulent as a result of both axial and crossflow 
viscous effects. In the numerical method axial motion is treated as if 
the flow were inviscid. For laminar flow, separation would begin on the 
leeward side of the body at an axial station upstream of the station at 
which turbulent separation begins. Also, the separation point would be 
closer to the windward side when the flow is laminar. From the pressure 
coefficient data for a circular cylinder (ref. 15), it was found that in 
turbulent flow at a Reynolds number of 6.7 x 10^, separation occurred at 
120°; in laminar flow at a Reynolds number of 1.85 x 10"*, separation 
occurred at 90°. Thus, the difference between the numerical and experi- 
mental separation positions is in qualitative agreement with the differ- 
ence between laminar and turbulent separation positions. To obtain more 
realistic theoretical predictions, it will, therefore, be necessary to 
generalize the equivalence principle to account for turbulent effects due 
to axial flow; it is felt that the principle, which was generalized for 
this program by including Newtonian viscosity in the crossflow equations, 
could be extended satisfactorily by the addition of eddy viscosity terms. 

4.4.4 Accuracy of the Ogive-Cylinder Calculations 

Based on the equivalence principle theory, an error analysis of the 
surface pressure results of figure 18 was made. The maximum slope of the 
ogive-cylinder body with respect to the free- stream flow direction was 
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28.92° which corresponds to a t = 0.552, and occurs on the windward side 
of the body at the nose. For this value of t, the hypersonic similarity 
parameter M t is 1.09, which satisfies the consistency condition of the 
theory. Since consistency is satisfied, the maximum absolute error in 
the surface static pressures, |AP|, that one should expect from this 
theory is about 19 percent of the free-stream dynamic pressure [see 
eq. (21)]. The maximum absolute error in surface static pressure, or 
maximum difference between numerical and experimental pressures, for each 
of the cross sections of figure 18 is presented in table 1. Table 1 
shows that the errors are very much less than the maximum error calcu- 
lated from the equivalence principle theory of Van Dyke. It is believed 
that the static pressure errors in the vicinity of the body's nose would 
be of the order of 19 percent and would decay rapidly as t decreased. 

For example, at the station 0.502 diameters, the parameter t is 0.481, 
which results in a predicted absolute maximum error of 11 percent of 
free-stream dynamic pressure. Therefore, the numerical results are con- 
sistent with the equivalence principle theory of Van Dyke. 

TABLE 1. MAXIMUM SURFACE PRESSURE ERROR IN 
VARIOUS CROSS-SECTIONAL PLANES 


Axial station, 
maximum cylinder diameter 


| AP | / q x 100, 
percent 


0.502 5 
2.990 1 
1.289 3 
3.970 5 
4.920 4 
5.830 4 
7.630 3 


4.4.5 Drag and Lift Coefficients 

Drag and lift coefficients were determined on the basis of numerical 
pressure data only; shear stress effects were not included. To compute 
these coefficients, the pressure coefficients on the ogive-cylinder surf- 
ace were numerically integrated to determine the lift and drag. 
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Comparison of the numerical and experimental coefficients, as shown in 
figure 20, was found to be very satisfactory. The coefficients of total 
drag from the two sources proved to be almost identical functions of 
distance from the nose of the ogive. Differences in the lift coefficient 
curves began to be significant at a distance of about 2.75 diameters 
from the nose of the ogive, increased to a maximum at about 4.5 diameters, 
and then decreased; the numerical lift coefficient was approximately 5.8 
percent higher than that measured at about 4.5 diameters from the ogive 
nose, and about 2 percent higher at 8.35 diameters from the nose. 

An overestimated lift coefficient is consistent with the fact that 
separation occurred earlier in the numerical flow field than in the 
experimental field. Lift coefficients are calculated from pressure 
forces on the surface of the body, as projected in a plane parallel to 
that of free-stream flow. Since the area of the appropriate projection 
is almost entirely derived from the body's long cylinderical surface, 
and separation takes place on that surface, separation that occurs too 
early increases the lift. 

On the other hand, the drag coefficient should be relatively inde- 
pendent of the location of separation. Since the ogive becomes cylindri- 
cal at a station 3.0 diameters from its nose, the projected surface area 
normal to the direction of free-stream flow is very small for the long 
cylindrical surface aft of this station, and the corresponding pressure 
forces contribute very little to the total drag. Quantitative agreement 
in the pressure drag coefficients is, therefore, physically reasonable. 

4.5 RESULTS OF FUSELAGE PROBLEM 

For the fuselage problem the crossflow Mach number was 0.288 and 
the crossflow Reynolds number was 1.15 x 10 /ft. Cross-section flows 
have been calculated in this problem in planes normal to the 7%° refer- 
ence line to a horizontal station 19.5 inches aft of the nose. At this 
station the curve bounding the body's surface represents an asymmetric 
fuselage configuration with a canopy. To reach this station the problem 
was run 2079 cycles (time steps) on the UNIVAC 1108 computer. This cal- 
culation required approximately 6 hours of computer time. 
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4.5.1 Qualitative Behavior of Numerical Flow Field 

From velocity vector plots of the cross-sectional flow field, the 
qualitative behavior of the fuselage flow can be examined. The sequence of 
events as one moves downstream along the 7%° reference line of the fuselage 
is indicated in figures 21 to 25. The vectors in these figures correspond 
to the particle velocities of the flow at each mesh point of the finite 
difference mesh; the tail of each vector corresponds to the mesh point. 
Initially, a bow shock forms at the nose and an expansion fan, caused by 
the interaction between the expanding body and the crossflow, appears on 
the leeward side. This is indicated in figure 21 at a station 1.135 inches 
from the fuselage nose. This flow pattern continues as the fuselage cross 
section grows until it reaches the canopy which induces a shock wave in the 
flow field (see fig. 23) . As the canopy radius reaches maximum (fig. 24) 
and starts to decrease, a rarefaction develops in the flow field above 
the canopy, leading to the formation of vortices (fig. 25) . 

4.5.2 Surface Pressure Comparison 

The quantitative behavior of the predicted flow about the fuselage 
was ascertained by comparing the numerical results to experimental mea- 
surements. The fuselage configuration of figure 6 was tested in the 
Ames 8- by 7-Foot Wind Tunnel. Static pressures were measured along the 
surface of the fuselage, and flow-field measurements with conical probes 
were made in cross-sectional planes normal to the norizontal reference 
line at a station 19.5 inches aft of the fuselage nose. The static- 
pressure instrumentation of the fuselage forebody surface is sketched in 
figure 26. The pressure taps were located in cross-sectional planes 
normal to the horizontal reference line. The point of intersection of 
the 7%° reference line with this cross-sectional plane defined the axis 
from which pressure tap locations were measured along the periphery of 
the cross section. Pressure taps were located -90°, -60° -30°, 0°, and 
+24° from this axis (see fig. 26). Thus, with respect to stations along 
the horizontal reference line, the instrumentation was located along the 
five planes indicated in figure 26. 
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Numerical and experimental distributions of surface pressure coeffi- 
cients are compared in figure 27 along the five planes of figure 26. It 
is seen that excellent agreement has been achieved along the -90°, -60°, 
and -30° planes. The numerical and experimental pressure coefficients 
agree along the 0° plane until station 18. Along the 24° plane there is 
a discrepancy at station 15. These discrepancies can be attributed to 
two sources: (a) a canopy shock wave-axial boundary-layer interaction, 

which occurs experimentally but not numerically, and (b) a slightly dif- 
ferent calculated canopy shock location which will be described in more 
detail in the discussion of the flow-field results. The canopy shock- 
axial boundary-layer interaction spreads the shock pressure rise over a 
greater distance than calculated. 

4.5.3 Accuracy of the fuselage Calculations 

As in the case of the ogive-cylinder problem, an error analysis was 
made of the surface static pressures. The maximum slope of the fuselage 
geometry with respect to the free-stream flow occurred on the windward 
side of the body at its nose (see fig. 8). Its value was 22°, which 
corresponds to t = 0.400 and M^t = 0.998. Therefore, the consistency 
condition of the theory is satisfied. According to equation (22), the 
maximum absolute error then becomes 5 percent of the free-stream dynamic 
pressure. On the basis of the results of figure 27, the maximum error, 

| AP|/q , in each of the planes was determined and the values are tabulated 
in table 2. It is seen from the table that the errors recorded are con- 
sistent with the maximum error calculated from the equivalence principle 
theory. 


TABLE 2. MAXIMUM SURFACE PRESSURE ERROR ALONG THE 
FIVE PLANES OF FIGURE 33 


:, deg 

|AP ]/q TO x 100 

90 

2.0 

60 

2.3 

30 

3.4 

0 

5.6 

24 

5.8 


percent 
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At first glance the relatively small errors of table 2 are surprising 
in lieu of the fact that hypersonic small disturbance theory is applied at 
supersonic speed. One would expect linearized supersonic theory to apply 
in this Mach number range. However, an analysis of the errors shows that 
linearized supersonic theory is only slightly more accurate than hyper- 
sonic small disturbance theory for the fuselage problem. The error in 
linearized supersonic theory is of the order t (ref. 2), which results in 
a maximum absolute error in the static pressure of 4% of free stream 
dynamic pressure. Therefore, for a relatively thick class of bodies 
(which includes most fuselage geometries), where the hypersonic similarity 
parameter is of order unity, both hypersonic small disturbance theory and 
linearized supersonic theory give satisfactory results. 

4.5.4 Contour Maps in the Flow Field of the Fuselage 

Contour plots of the local pitot pressure ratio (P /P^) , local 

sideslip angle (P) , local Mach number (M) , local angle of attack (a g ) , 

and local total pressure ratio (P /P ) were generated from the numerical 

s s°° 

results and compared to corresponding contour plots from experimental 
data. The comparisons were made in a cross-sectional plane normal to the 
horizontal reference line at a station 19.5 inches from the fuselage apex. 
See figures 28 through 32. The above parameters are defined in Appendix E. 

The local pitot pressure, P , is effectively the stagnation pressure 
measured by a probe whose axis is parallel to the local flow direction. 

As in the case of the surface pressures the pitot pressure is measured 
directly; therefore, it provides a reliable measurement for comparison. 

The numerical pitot pressure was calculated from the Rankine-Hugoniot 
equations for a normal shock, where upstream of the shock, the local Mach 
number, and stagnation pressure were assumed to exist. The numerical and 
experimental pitot pressure ratios (P /P ) of figure 28 indicate good 

p Soo 

agreement in the lower quadrant of the flow field, where the two sets of 
P /P contours are nearly coincident in the vicinity of the body. This 

P S CD 


34 



quantitative agreement is in accord with the agreement obtained between 
surface static pressures on the -30° plane, which intersects the body in 
this region. In the upper quadrant of the flow field there is a 
discrepancy between the experimental and numerical contours. The numeri- 
cal results indicate that the canopy shock intersects the fuselage surface 
above that point indicated by the experimental data. The increased down- 
ward distance of travel by the canopy shock in the experimental case 
accounts for the discrepancy in contours in the upper quadrant of the flow 
field. Although axial-boundary layer effects and discretization errors 
influence the canopy shock location, it is believed that the equivalence 
principle approximation of neglecting the velocity perturbation along the 
7%° reference line is the primary cause of this discrepancy. The velocity 
perturbation is effectively zero upstream of the body apex, hence, the 
windward oblique shock emanating from the apex was calculated correctly. 
However, the velocity perturbation is positive upstream of the canopy, 
and since this quantity is neglected in the method, the calculated canopy 
shock must be, and is, weaker than experimentally measured. 

The local sideslip angle, g, is measured in the cross-sectional 
plane normal to the horizontal reference line and is the flow inclination 
in the x',z' plane (see fig. 9). As in the case of the pitot pressures, 
the comparison of the experimental and numerical sideslip angle contours 
shown in figure 29 indicates that the predicted flow field in the lower 
quadrant is nearly correct. Furthermore, the predicted location of the 
canopy shock and the predicted flow-field contours in the upper quadrant 
deviate from experiment as in the pitot pressure comparisons. 

Contours of constant Mach number and constant local angle of attack 
with respect to the horizontal reference line (see fig. 8) are presented 
in figures 30 and 31, respectively. It is seen from these figures that 
the numerical Mach number and angle-of-attack contours do not match the 
experimental contours even in the lower quadrant of the flow field. This 
discrepancy is not in accord with the surface static pressure comparison, 
pitot pressure comparison, and sideslip angle comparison discussed 
previously. 
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The reasons for this disagreement can be found in discretization 
errors in the numerical method, errors introduced by the equivalence 
principle assumption, and errors inherent in the experiment and data 
reduction process. The error in Mach number due to the equivalence prin- 
ciple has been determined previously [see eq. (28)]. Equation (28), 
evaluated with t = 0.40 and M = 2.5, results in a maximum absolute error of 

CO 

about 30 percent of the free-stream Mach number. This corresponds to an 
absolute Mach number error of about 0.75. Consider the experimental 
contour of figure 30 for a constant Mach number of 2.50. It is seen from 
figure 30 that this contour coincides with the numerical contour for a 
constant Mach number of 2.40, which results in an absolute Mach number 
difference of 0.10. Thus, the observed absolute difference in Mach num- 
ber between the numerical and experimental results is within the 
maximum error one would expect from the equivalence principle alone. 
Therefore, one must conclude that errors introduced by the equivalence 
principle assumption account for the major part of these discrepancies. 

The calculated local total pressure recovery contours (P /P ) are 

1 s S co 

shown in figure 32. The predicted total pressure recovery varies from 
1.0 to 0.88 in the flow field, with recoveries near 1.0 throughout most 
of the flow region. This parameter is the most difficult to calculate 
accurately and to determine experimentally. Although the experimental 
data for this fuselage configuration are still preliminary, an error 
analysis indicates that the values of P /P can be determined only with 
± 0.03. Por a total pressure ratio bandwidth from 0.88 to 1.0, this 
error is too large to yield any significant contour data. Consequently, 
no experimental data are shown. Future refinement of the set of conical 
probe calibration data will hopefully reduce this error to a more 
meaningful level. 
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4.5.5 Lift and Drag Coefficients 

Numerical lift and drag coefficients for the fuselage problem are 
compared to corresponding lift and drag coefficients determined from 
inviscid, linearized supersonic theory for flow about an axisymmetric 
body at angle of attack (ref. 16) . Linearized supersonic theory was used 
because there were not enough pressure data to determine the lift and 
drag experimentally. The lift and drag coefficients from the numerical 
method were determined by numerical integration of the surface pressure 
coefficients over the surface and are based on a cross-sectional area 
corresponding to the maximum equivalent radius of the fuselage, 
a max = ^ inches. As in the case of the ogive-cylinder problem, shear 
stress effects were not included. It is seen from figure 33 that the 
numerical lift coefficient distribution nearly corresponds to that from 
linearized supersonic theory. In view of the fact that the maximum 
absolute static pressure errors are about the same for both linearized 
supersonic theory and hypersonic small distrubance theory, agreement in 
the lift coefficient is expected. On the other hand, the numerical drag 
coefficient distribution shows greater drag than the distribution from 
linearized supersonic theory. This discrepancy may be attributed to the 
fact that the coefficients determined from linearized supersonic theory 
have been further specialized to the case of very slender bodies (ref. 16) , 
which is not true for this geometry. However, this comparison does 
indicate that the lift and drag coefficients calculated by the numerical 
method of this paper, if not correct, are at least of the right order. 
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5.0 CONCLUSIONS AND RECOMMENDATIONS 


The principal conclusion emerging from the work reported here is 
that a promising new method has been developed for calculating flow fields 
about slender bodies of arbitrary geometry in the supersonic and hyper- 
sonic flight regimes. The present version of the numerical method can 
adequately predict the surface pressures and flow field characteristics 
outside the boundary layer for bodies where the hypersonic similarity 
parameter M^t is unity or greater and for free stream Mach numbers 2 and 
above. For cases where boundary layer separation and the subsequent 
formation of spiral vortex sheets occur, the numerical method predicts too 
early a separation which results in small inaccuracies in predicted flow 
field quantities downstream of the region of separation. 

The basis for the above conclusion lies mainly in the specific 
results obtained for the ogive -cylinder and the fuselage configurations. 
For the ogive-cylinder problem a spiral vortex pair was computed on the 
leeward side of this body. Although axial effects caused some differ- 
ences between numerical and experimental separation regions and vortex 
center positions, in general, the static pressure distributions on the 
body's surface, the lift, and pressure drag of the body agreed with the 
corresponding experimental values; the surface pressures differed by no 
more than 5 percent of free-stream dynamic pressure; the lift coefficients 
differed by no more than 6 percent; and the drag coefficients by no 
more than 2 percent. From velocity vector plots of the cross-sectional 
flow field of the fuselage configuration, the structure of the flow field 
seemed correct, at least qualitatively. Furthermore, comparisons of the 
numerical static pressure distributions on the fuselage surface were in 
quantitative agreement with experimentally determined distributions; the 
maximum error in the static pressure was no more than 6 percent of free 
stream dynamic pressure. Contour plots of pitot pressures were of the 
same general shape as the experimental contours and some of these curves 
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coincided. Finally, the pressure drag and lift distributions seemed 
correct when compared to corresponding drag and lift coefficients derived 
from inviscid, linearized supersonic flow theory. 

The justification for the part of the above conclusion which per- 
tains to the range of applicability of the method, i.e., M^t s 1 and 
M S 2, lies mainly in the fact that for both problems (M t » 1. M = 1.98, 

CD CO ? CO 

and M =2.5) the deviations predicted from experimental results were 

CO 

contained well within the predicted error bounds of the equivalence prin- 
ciple. Hence, as the Mach number increases, the Van Dyke theory tells 

us that the predicted results will become more accurate. However, to 
strengthen this conclusion, additional comparisons should be made between 
numerical and experimental data at other flight conditions and for differ- 
ent geometries within the above range of applicability. 

To realize the full potential of this numerical method, we believe 
that further research is justified. Three areas of additional research 
seem promising. First, the method should be revised and include the 
velocity perturbation on the principal axis along which the equivalence 
analog is made. The inclusion of the axial velocity perturbation in the 
finite difference equations of motion will improve the accuracy of the 
method in its present range of applicability and will extend the method 
to cases where M t « 1 and for Mach numbers above transonic. Hence, a 

CO 

unified method, which applies throughout the supersonic and hypersonic 
flight regimes, will be achieved. 

The second area where additional research seems justified concerns 
the relationship between viscous axial and crossflow effects. In the 
ogive cylinder problem, too early a separation was calculated considering 
only viscous crossflow terms which were assumed laminar. The discrepancy 
in the separation point location is attributed to both axial viscous 
effects and turbulence. However, from the results of this problem it was 
not possible to determine which was the more important effect. There- 
fore, to determine the relationship between viscous axial and crossflow 
effects it is necessary to solve an additional problem (both numerically 
and experimentally) at hypersonic flight conditions and in the laminar 
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flow regime. At hypersonic flight conditions, errors in the equivalence 
principle are minimized; thus, deviations between numerical and experi- 
mental results can be attributed entirely to the neglect of axial viscous 
effects. If these deviations are large, some of the more important axial 
viscous terms can be retained in the finite difference equations of 
motion. 

Finally, turbulence effects can be included in the finite difference 
equations of motion by the addition of eddy terms, i.e., eddy viscosity 
and thermal conductivity. These eddy terms would be based on empirical 
constants which would be determined from comparisons of numerical and 
experimental results. 

Along the lines just sketched, we believe that an efficient, 
reliable capability can be developed for predicting viscous flow field 
characteristics for nonsimple geometries throughout the supersonic and 
hypersonic flight regimes. 
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APPENDIX A 


THE EQUIVALENCE PRINCIPLE 

In this section the equivalence principle for steady, inviscid 
three-dimensional flow is derived. An Eulerian coordinate system fixed 
with respect to the body and having its z-axis along the axis of the body 
is used throughout this section (see figure 2) . 


The equations of motion for steady, inviscid three-dimensional flow 
about a body are: 

Continuity 


(Pu) + (Pv) + ^ [P (U^Cos a 4- w)] - 0 


(Al) 


Momentum 


du du . s du 1 dP 

u ^ h v ^ f- (w + U Cos <a) = - — 

dx dy v 00 ' dz n 


p dx 


(A2a) 


dv dv . dv 1 dP 

dx dy K °° ' dz P dy 


(A2b) 


dw , dw , . , _ . dw 

U d^ +V d^ + < w + U =° Cos «> di 


1 a?. 

P dz 


(A2c) 


First Law 


u|^ + v|^ + (w + U Cos O') = 

dy 00 dz 


9x 


(A3) 


- P 


r d /] 

l\ , d / 1 \ 

L u dx (f 

>) V dy (p) 


(?) + <* + D » Cos a) h (?)] 


41 


Boundary Conditions at Body 


(U^Cos a + w) 


SL + v 4 + u 4 = 0 

z y x 


(A4a) 


Boundary Conditions at Shock 


p U (Cos an + Sin a n ) = p [(U Cos a + w ) n +vn + u n ] (A4b) 

“ z y s L 00 s z s y s x J 


2 2 
+ P m U m (Cos a n + Sin an) = P 

' co ' co co x 2 y g 


+ p f(U Cos a + w ) n +vn + u n T 
r s L “ s z sy sxj 


(A4c) 


P 

„ , + -r U" (Cos an +Sinan ) 2 = — Y - — 

Y “ 1 Pa. 2 z y Y - 1 P s 


jy 1 2 


+ i f(U Cos a + w ) n + v n + u nl 
2 L v “ s z s y s xJ 


(A4d) 


where E is the internal energy per unit mass, P is the pressure, p is the 

density, U ro is the free stream speed, w is perturbation velocity in the 

z-direction, v is the velocity in the y-direction, u is the velocity in the 

x-direction, and 4,4, and 4 are the direction cosines of the body normal 
x y z 

in the x, y, and z directions respectively. The subscript s refers to 

properties downstream of the shock, and n , n , n are the direction cosines 

x y z 

of the shock normal in the x, y, and z directions respectively. 

Boundary Conditions Far From Body 


z = - o= : w = U Cos a, v = USina, u = 0 

CO 3 CO J 


(A5) 


Introducing the transformation equations 


z = U^Cos at 

y = y 


(A6) 
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and the slender body assumption 


u » w 

CO 


into equations (Al) through (A6) yields 
Continuity 


Momentum 


t£ + b (pu) + h (pv) = 0 


3w , 3w , 3w 

rr + u r 1- v r— 

ot ax ay 


Ik/ 1 \ 

p Bt \ u Cos a > 


First Law 


3v , Bv . Bv 1 Bp 

tt— • + u ~ — + V •r — = - — tr*- 

3t 3x By p By 


£ii + u ki. + v 3u = _i<£ 

3t 3x By p Bx 


BE , BE BE 
Bt + u Bx +v 3y = ” P 


h (?) + u h (p ) + v 1? (?) 


Boundary Conditions at Body 


Cos a U & + u £ + v j £ 

co z x y 


0 


Shock Conditions 


p^U^C- Cos a Sin p + Sin a Cos p) 


= p (v n 
s s y 


+ u n 
s x 


U Cos a Sin P) 
00 y 


(A7) 

(A8a) 

(A8b) 

(A8c) 

(A9) 

(A10) 

(Alla) 
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(Allb) 


11111 


P oo + Poo^oo^" CoS a Sin P + Sin a Cos P) 2 = P s 

2 

+ p(vn + u n - U Cos o' Sin P) 
s sy sx °° ' 

y P co U *(- Cos a Sin P + Sin a Cos P ) 2 Y P S 

V - 1 + 2 = Y - 1 pj 

+ ( v s n y + u s n x - U^Cos a Sin P) 2 (Allc) 


where: P is the angle made by the shock surface with respect to the z-axis 

in the (y,z) plane. Equations (A7) , (A8b) , (A8c) , and (A9) are independent 
of w and, thus, represent the equations for time-dependent motion in the 
(x,y) plane. Therefore, the steady 3-D equations of motion have been trans- 
formed to a time-dependent 2-D set of equations. The 3-D shock boundary 
conditions also reduce to a nonsteady shock in the (x,y) plane moving with 
velocity proportional U^Cos a Sin P. In summary, the necessary and 
sufficient conditions for the equivalence principle to be valid are: 

a. The z-component of the local velocity vector must be 
approximately equal to the z-component of the free stream 
velocity vector, and 

b. the time -dependent solution in the (x,y) plane must satisfy 
the three-dimensional boundary condition equation (A10) at 
the body surface. 

The calculational procedure for inviscid flow is then as follows: 

First, equations (A7) , (A8b) , (A8c) , (A9) , and (A10) are solved for v,u,P, 
and p. The perturbation velocity w may be subsequently obtained from the 
Bernoulli equation, which is 


Y 


Y 


1 


(w + U^Cos a) 2 + v 2 + u 2 
P + 2 


constant 


(A12) 


Therefore, the equivalence between a steady- three-dimensional flow and a 
time-dependent two-dimensional flow is proven. 
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If the conditions (a), (b) are satisfied, viscous effects can be 
included in the time -dependent calculations in the (x,y) plane without 
violating the equivalence principle assumptions. A no-slip boundary 
condition can be applied at the surface of the cross section in the 
(x,y) plane in conjunction with the equivalence principle boundary 
condition equation (A10) . 
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APPENDIX B 

CONSERVATION OF TOTAL ENERGY 


In this section the self-consistency property of form that the 
finite difference equations (8) to (13) possess is demonstrated. This 
demonstration proceeds in three steps. First, the finite difference 
kinetic energy equation is derived for the momentum zone a^ 1 / 2 ’ 

3 4+1/2’ Then t ^ le finite difference First Law equation is derived for 
this momentum zone. Finally, these relations are added to determine 
the finite difference equation of total energy for the momentum zone 

a j&-l’ 3 4- 1/2 " 

The kinetic energy equation for momentum zone a. ^ ]_/2* 3 4+1/2 

can be derived from momentum equation (12) which is centered at the 
1 / o 

time t . Equation (12) can be written in terms of velocities by 

employing the forward extrapolation relation, equation (13), i.e., 

(u n+1/2 - u 11 ' 372 ) / \ 

\ U 4 ) .. /pH- 1/2 . pH- 1/2 \ 


2 m 4 


At 


4 - 1/2 


4+1/2 J 


(Bl) 


Multiplication of equation (Bl) by U 


n-1/2 

4 


yields the kinetic energy 


equation for momentum zone a^ -^/2 > 3 4+1/2’ 


At 


U 


n+1/2 n-1/2 


4 


4 


= U 


n-1/2 / p n-l/2 

\ P 4-1/2 


n- 1 / 2 n-3/2 


n- 1/2^ 
^4+1/2/ 


(B2) 


The First Law equation for momentum zone a^ 3 4+1/2 is 

derived as follows: the internal energy for this zone equals half 

the sum of the internal energies of thermodynamic zones a^ a ^ and 
a 4., a 4+l ( see equations (9) and (10)) . This division of internal energy 
comes directly from the relationship between thermodynamic and momentum 
zones specified when the physical model of figure 6 was postulated. The 
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-P^r term for momentum zone a SL+l/2 comes from half the - 

for thermodynamic zones , a^ and a^, a The First Law equation 

for momentum zone a ^_^/2’ a JlH-l/2 ^ ecomes: 


2 W- 1/2 + E m/ 2 ) - l( E £/ 2 + C 1 / 2 ) 


1 ( P £-l/2 + P .fc- 1 / 2 ) / TT n-l/2 T7 n- 


- u n " 1/2 ) 
U >1 / 


^ p a ” ^ + p n ^ 

1 1 + 1/2 m / 2 / / n- 1/2 T n- 1 / 2 \ 

2 \ U X+1 " / 


The finite difference equation for conservation of total energy 
for momentum zone a^^, a f+l/2 resu - Lts from the addition of equations 
(B2) and (B3) . 


( H ! ' H £ ) / n- 1/2 n-1/2 \ 

— T r>i /2 " ™ 1 - 1 / 2 ) 


where: 


un-1/2 un+1/2 

u n _ / E n + E n \ I + Jl 1 

~ rA-1/2 >1/2/ 2 2 

L/2 ( f Il/2 + ^j/z) ("I 17 " + 

2 


n-1/2 _ li-1/2 
> 1/2 

(P n 

n-1/2 = r 1+1/2 
>1/2 


P ll / 2 ) (' 


T n- 1/2 . ,,n-l/2 


Equation (B4) , derived from finite difference analogs of the 

First Law and momentum equations, is a reasonable finite difference 

expression for total energy conservation. The rate of work done on 

n 

the surface, having the Lagrangian coordinate & at t: *- me t , 

is the product of the time-averaged pressure in the thermodynamic 
zone a^ a^ and the space-average velocity between surfaces a ^ ^ 
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and a^. The internal energy of the momentum zone a^ y_/2* a M-l/2 

at time t n , is half the internal energy of the two thermodynamic zones 

which contain it. Finally, the kinetic energy of momentum zone 

a . 1 a. at time t n , is the product of the velocities at the 

Xi-i/z xrhi/z n-1/2 n+1/2 

surface a^, at times t and t . It is believed that this self- 

consistency of form property, which the differential equations possess, 

and which the finite equations preserve, is the primary reason for their 

success in numerical calculation of one-dimensional time-dependent flow 

fields . 
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APPENDIX C 


AFTON 2 PE FINITE DIFFERENCE EQUATIONS AT INTERIOR MESH POINTS 


The AFTON 2PE finite difference equations are presented for a 
generalized coordinate system. The field of motion in AFTON 2P is actually 
covered with two closely related finite difference meshes -- one for the 
calculation of thermodynamic variables such as stress, i.e., quadri- 
lateral zones, and the other for the calculation of kinematic variables 
like momentum, i.e., momentum zones (refs. 11 and 12). Figure 4 illustrates 

the two types of meshes in two space dimensions. The continuity and first 
law equations are applied to calculate properties on a quadrilateral zone 
while the equations of momentum and total energy conservation are used 
to calculate properties on a momentum zone. In the forthcoming analysis 
quadrilateral 1, 2, 3, 4 and momentum zone a, b, c, d are considered 
(see figure 4). The finite difference equations in the generalized 
coordinate system are as follows: 

Definitions : 


R = R(x,y) Cartesian coordinate position relative to the 

moving frame 

U = U(u,v) material velocity relative to the laboratory 

frame 


E 

i 


H 


internal energy 

unit vector in the x-direction (i.e., the plane of 
flow, and normal to the free stream flow relative 
to the laboratory frame) 

total energy 


W 


rate of work 


S = S/x,y) mesh point velocity relative to the laboratory 

frame 


P 


material density 
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m 

V 

M 


fJ- 

t 

At 

y 


xy 


O' 


it 


( )° 
( )' 
( )" 

( ) 

( ) 


unit vector normal to (x,y) plane 

mass 

volume 

momentum 

viscosity 

time 

times tep, t = t' - t° 
ratio of specific heats 
fluid pressure 
viscous stresses 



denotes property at initial time, t° 

denotes property at final time, t' 

denotes property a half timestep before initial 
time t 

denotes property a half timestep before final 
time t ' 

denotes timestep, a half timestep after final 
time t' 


Equations 


R' = R'(jt', cross sectional geometry) 
R = % (R ' + R°) 

S = (R' - R°)/At 


(Cl) 

(C2) 

(C3) 
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U = U(u, v) 


—12 

-32 

-12 

—12 

—32 

-32 


(pWA) 

(pWA) 


12 

32 
V' 

V 
1 

m = m 


d. 

IX 


(% " ^ x iS 

(R 3 - R 2 ) x k 

+ u 2 ) 

%CS X + s 2 ) 

%(U 3 + U 2 ) 

%(S 3 + s 2 ) 

P 12®12 " - 12 ) 

P 32®32 " — 32^ 

V'(t', cross-sectional geometry) 

(V ' + V°) 

° +At[-(pWA ) 12 + (pWA) 32 + ( P WA ) 43 - (pWA) 4] J 

1 /17 1 
m /V 

X i( y i -1 " y i+ 2 ) + X i-l( y i +2 " y i) + X i +2 ( y i - y i-l) 
where: i=l,2,3,4 and, for example, when i=l , i-l=4 


-12 

-32 


(C4) 

(C5) 

(C 6 ) 

(C7) 

(C8) 

(C9) 

(CIO) 

(Cll) 

(C12) 

(C13) 

(014) 

(C15) 

(C16) 

(C17) 


iy 

A. 

IX 


it 


0U 

3x 


-d. 

IX 


% [ (y i-l ' y i+2 )/d ix + (y i-l ' y i+l )/d i+l x 

+ (y i+2 - y i-l^ /d i + 2 x] 
" X i+2 )/d iy + (X 1-1 ’ X i-l )/d i+l x 

+ < x i+ 2 ' x i-l )/d i+2 y] 


A lx U l + A 2x U 2 + A 3 x U 3 + A 4x U 4 


(C18) 

(C19) 

(C20) 

(C21) 


—■ = A. u + A„ u + A„ u„ + A. u 
3y ly 1 2y 2 3y 3 4y 4 


3v _ 


A. v, + A„ v„ + A„ v n + A. v. 


9x lx 1 2x 2 3x 3 4x 4 


(C22) 

(C23) 


= A, v. + A„ v + A„ v„ + A. v. 
3y ly 1 2y 2 3y 3 4y 4 


(024) 
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•;-*[!!?-!£] (c25) 

<c26) 

^ --*[%+ S] <c27) 

P = (V - 1) P 1 E° ( c28 ) 

\ 2 = h(\ ■ R 2 ) 4 (C29) 

1^2 = < 7 ° + p - (a") 1 - ( c3 °) 

(pWAE) 12 = ( P WA) 12 E 12 ( C31 > 

E t = - (pWAE) 12 + (pWAE) 32 + (pWAE) 43 - (pWAE) 41 (C32) 

i . [«°e° - (^-Dj + £ 13 '2 2 + £ 24 '2 3 + Eaa-U* < c33 > 

- E t )AtJ /m 1 

P = (Y - 1) p 1 E ( c34 ) 

% 2 = %(p° - p) • 4 2 (c35) 

E = E - [At(| 42 -U 1 + f 13 -U 2 + F 24 -U 3 + l 31 'V] M 1 ( c36 > 

P = (Y - 1) P 1 E ( C37 ) 

2=^0° + tr 1 ) • ^ ( C38 ) 

E 1 = [m°E° - + F 13 -U 2 + + P 31 -ty&t] /m 1 (C39) 

m'*' = + tn^ +4 + m d) (C40) 
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M fc = %[( P WA) 42 + (pWA) 12 + (pWA) 54 + ( P WA) 61 ] • 

+ %^(pWA) 12 + (pWA) g9 + (pWA) 61 + (pWA) 7g J • %(U X + Ug) 

+ %[(pWA ) 56 + (pWA ) 41 + (pWA) 6? + ( P WA) 18 ] • %(Ug +U X ) 

+ %j^(pWA ) 41 + (pWA) 82 + (pWA) lg + ( p WA ) 29 j • \ (_U^ + U 2 ) 

M 1 =M° + (^2 + Fg 4 + £*> +F 2g +MP At 


= (2M 1 /mp) - U 


1 . _1 . 1 „1 , 


H 1 = inn* (U x * U x ) + %(< E a + “b E b + m c E c + m d E d } 

4d = a ‘ ( ~a • -V X £ 

W = - [%( U x + U,) • (4 b _ + F^ c ) + ^ + Ug) ■ (F c£ 4- 

+ %<£, + U 2 ) • (F^ + Ig a ) +h(V 1 + V • (— aa + — 

Hg = %[(pWAE) 43 + (pWAE) 12 + 0?WAE) 59 + #>WAE) 6 J 

+ \ |^(pWAE) 12 + ( P WAE) 89 + (pWAE) 61 + (pWAE)^] 

+ %^( P WAE) 56 + (pWAE) 41 + (pWAE) 6? + (pWAE) lg ] 

+ %j^(pWAE) 41 + ( P WAE) 32 + (pWAE) lg + (pWAE) 29 J 

Hre = Te [(P WA ) 43 + <P WA) 12 + ( P WA) 54 + ( P WA) 6l] + 

+16 [(pWA) 12 + ( P WA) g9 + (pWA) 61 + ( p WA) 7g ] + U 8 ) 2 

+16 [<P WA) 56 + (pWA) 41 + (pWA) 67 + ( P WA) ia] ( ^6 + ^ 

^6 [(pWA) 41 + (pWA) 32 + (pWA) lg + (pWA) 29 ] + U 2 ) 2 


(C41) 


(C42) 

(C43) 

(C44) 

(C45) 

(C46) 

(C47) 


(C48) 
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APPENDIX D 

PARAMETRIC REPRESENTATION OF FUSELAGE CROSS SECTIONS 


In this section the parameters which describe the fuselage cross sections 
are specified as functions of distance along the 7%° reference line (see 
f igure 8) and the 7%° reference line is prescribed as a function of horizontal 
station. The cross sections taken were perpendicular to the 7%° reference 
line (i.e. z-axis, figure 8). The parameters of each cross section was deter- 
mined from geometrical data which were taken from the design blueprint of the 
fuselage geometry. 


A schematic representation of the fuselage cross section normal and the 
7%° reference line is presented in figure 9. In figure 9 the parameters which 
describe the cross section are indicated. These parameters are defined as 
follows : 

(A^(z), B^(z)), coordinates of the center of the smaller circular 
arc in the fuselage periphery 

R^(z), radius of the smaller circle of the fuselage periphery 

(A^Cz), B^Cz)), coordinates of the center of the larger circular arc 
in the fuselage periphery 


R^Cz), radius of the larger arc 

(A,_(z), B^(z)), coordinates of the center of the arc of the canopy 

Rg(z), canopy radius 

In the case being considered, B^Cz) = A^(2) = A^(2) = 0, R^Cz)- R^(z) = i^.(z) 
(Since the shape supplied has flat sides and bottom), and R 2 (z)- B^(z) ^ R^(z). 
Furthermore, B,-(z) = R 2 (z) and R^(z) equal zero for z less than z^ 5 where Z2(=10.88") 
defines the distance along the 7% reference line at which the canopy first 
appears. The fuselage cross-section is circular to z^ = 9.44 inches and 
asymmetric between z^ and It proved convenient to define these parameters 

by another set of functions from a distance of z^ = 15.05 inches onward. The 
equations for the above parameters considered units of inches and are as 
follows : 
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(i) Radius of smaller circle of fuselage periphery 


12 R 


For 0 < z < 30.68, R, (z) = V a T 

1 Pi J 


1 J 
z 


where a 1 1 = (1.7873231 x 10 -1 ) 


l -? 

= (7.0559487 x 10 ) 

R 1 -2 

a 3 = (-3.2209788 x 10 ) 

R 1 -2 

a 4 = (6.6039801 x 10 ) 

R 1 L 

a 5 = (-7.5565592 x 10 - ^) 

R 1 -5 

a g = (5.1130544 x 10 


^ r 

a ? = (-2.1016895 x 10 °) 


R. 


- 8 , 


a g = (5.1489881 x 10 u ) 

ag 1 = (-6.6633139 x 10 -10 ) 
R 

1 10 
R ] 

a n = ( 4 ‘ 7027051 x io” ) 

a 12 = (“ 4 * 0954467 x IO " 16 ) 


1 -12 
a = (2.0697609 x 10 ) 


(ii) Radius of larger arc 


12 R 

For 0 < z < z 3 , R 2 (z) = a j 

J=1 


2 J 
z 


where R 2 

a. 



(3.4725590 x 10 -1 ) 
(-1.0362533 x 10 _1 ) 
(3.6798959 x IO -2 ) 
(-7.4940700 x IO -3 ) 
(9.1822566 x IO -4 ) 
(-7.2186739 x IO" 5 ) 


2 ft 

a y = (3.7708718 x 10 ) 

R 2 -7 

a g = (-1.3236585 x 10 ') 

R 2 -9 

a g = (3.0861620 x 10 ) 

R 2 -11 
a 10 = (“ 4 - 5831105 x 10 ) 

R 2 -13 

a n = (3.9237683 x 10 

R 2 -15) 

a 12 = O 1 - 4741591 x 10 J 
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(-4.1445421 x 10“ 6 ) 
(1.1565654 x 10" 7 ) 
(-1.9920283 x lO -9 ) 
(1.9332318 x 10" 11 ) 
(-8.0938014 x 10 -14 ) 



APPENDIX E 


DATA REDUCTION COMPUTER CODE 

A data reduction computer code has been developed to reduce the 
flow field data of Problems 211.0 and 212.0 to meaningful parameters. 

The function of the code is to compute the Mach number, flow angularity, 
total pressure recovery, and pressure coefficient in cross section planes 
normal to the body axis. These parameters are to be calculated from the 
more fundamental AFTON 2PE numerical solution data (i.e., pressure, 
density, velocity, etc.). 

The data reduction code, called "AMSD", performs three functions. 

First, it selects and reads a dump (i.e., the flow field properties at a 
particular axial station) from the AFTON 2PE dump tape. Then it calcu- 
lates the specific internal energy (E) , density (p), pressure (p) , and 
the velocity component (w) in the axial direction at each of the mesh 
points. The specific internal energy, density, and pressure are zone 
centered in AFTON 2PE, so the mass and internal energy of the momentum 
zones must first be evaluated in order to determine the specific internal 
energy and density at a mesh point (ref. 10) . The pressure is then determined 
from the perfect gas law. Based on these parameters and the velocity com- 
ponents u and v in the plane of calculation, the velocity component w in 
the axial direction can be evaluated from the Bernoulli equation as follows: 


w 


U cos a + 

CO 


K 




v 3 + u 3 
2 


U cos a 

CO 


(El) 


where 

U free stream velocity 

CO 

a , angle of attack of free stream flow with respect to body axis 

y , ratio of specific beats [y = (C^/Cv) ] 

P , pressure 

p , density 
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K = J^l + Sin 2 O' MwJ ( E2 ) 

where: 

a», free stream sound speed 
Mo, free stream Mach number 

Equation (El) is valid for a weak shock and where viscous effects are small. 

In the ogive-cylinder problem at an axial station of .502 cylinder diameters 
and at a point on the dividing streamline 1.36 cylinder radii upstream of the 
cylinder's center, the ratio of the local stagnation pressure (P) to the free 
stream stagnation pressure (P ) is .92. For the fuselage problem total 
pressure recovery Eg/E sm varies from .90 to 1.00. It is shown below that the 
entropy change is directly related to the total pressure recovery; hence, the 
entropy change is small in both problems and the weak shock approximation is 
valid. 

Based on the properties just cited, the required parameters are determined by 
the AMSD computer code at each of the mesh points, from the following relations: 


M r 2 u. 2 2 *1^ 

M = Lw + u + v J 

c 

c = [y(y-D e] % 


Ps 
Ps co 


= exp 


(S-Sqq ) . Y , 

f C Vr 


where : 


S-S<a 

c 


log 


<£=>'<& ^ 




hoo = Eoo + ~~ 

Pco 

and finally the pressure coefficient is evaluated from 

c - 

^ %Pos Uoo 


(E3) 

(E4) 

(E5) 

(E6) 

(E7) 

(E8) 

(E9) 
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where : 


C, local sound speed 

C , local pressure coefficient 
P 

h, local enthalpy 
S, entropy 
M, local Mach number 

Subscript infinity near any property refers to free stream conditions. 


The local angle of attack of the flow, o^, and the local angle of side - 
slip of the flow, 0, are measured in the coordinate system consisting of the 
horizontal reference line, the normal line to it in the plane of symmetry of 
the body, and the coordinate normal to both lines, i.e., x', y' s z' indicated 
in figure 9. On this basis and 3 are defined as follows: 


O' 

e 


3 


-1 \ 
tan ■ y ) 

(E10) 

z 

1 u 
.an’ 1 

(Ell) 


where: 


U* = u 

Uy = v Cos 9 + w Sin 9 

U 2 = -v Sin 0 + w Cos 9 

9 » tan^B^Cz' ) 

(f 

B^Cz') = local slope of central axis with respect to 
horizontal reference line 


The pitot pressure, Pp, corresponds to the stagnation pressure measured 
by a probe whose axis is placed parallel to the local flow direction. An 
approximation of the actual pitot pressure can be obtained by causing the 
local velocity vector to go through a normal shock. The relation is as 
follows: y 
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APPLICATION OF THE 

EQUIVALENCE PRINCIPLE TO AN AX I SYMMETRIC BODY 


SHOCK CROSS SECTION 



EQUIVALENT 

THREE DIMENSIONAL STEADY PROBLEM TWO DIMENSIONAL UNSTEADY PROBLEM 

FIGURE 2 



FLOW FIELD ABOUT AN AXI SYMMETRIC BODY 
AT ANGLE OF ATTACK 



FIGURE 3 



FLOW REGIONS ABOUT A BODY MOVING AT 
SUPERSONIC VELOCITY 



FIGURE 4 


ONE-DIMENSIONAL, TIME -DEPENDENT MOTION 
IN LAGRANGIAN COORDINATES 


PARTICLE TRAJECTORIES IN A 
LAGRANGIAN COORDINATE SYSTEM 



THERMODYNAMIC AND MOMENTUM ZONES FOR 
ONE-DIMENSIONAL TIME-DEPENDENT MOTION 
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FIGURE 5 




SCHEMATIC DIAGRAM OF A FINITE DIFFERENCE MESH 
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FIGURE 6 



OGIVE-CYLINDER BODY GEOMETRY 





SCHEMATIC OF FUSELAGE CONFIGURATION 



FIGURE 8 
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DEFINITION OF FUSELAGE CROSS-SECTIONAL PARAMETERS 
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FIGURE 9 



FINITE DIFFERENCE MESH FOR OGIVE- CYLINDER 

RADIUS O.OI ft 


CROSS FLOW 



FIGURE 10 



FINITE DIFFERENCE MESH FOR OGIVE- CYLINDER 

RADIUS 0.046875 ft 



FIGURE 



FINITE DIFFERENCE MESH AT STA 7.0 



FIGURE 12 



FINITE DIFFERENCE MESH AT STA 25.0 


CROSS FLOW 



FIGURE 13 
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VELOCITY VECTOR PLOT OF FLOW FIELD 
0.502 MAX BODY DIAMETERS FROM NOSE 
RADIUS 0.0145 ft 



FIGURE 14 


VELOCITY VECTOR PLOT OF FLOW FIELD 
2.99 MAX BODY DIAMETERS FROM NOSE 
RADIUS 0.046875 ft 



FIGURE 15 



VELOCITY VECTOR PLOT OF FLOW FIELD 
4.92 MAX BODY DIAMETERS FROM NOSE 
RADIUS 0.046875ft 



FIGURE 16 



VELOCITY VECTOR PLOT OF FLOW FIELD 
8.35 MAX BODY DIAMETERS FROM NOSE 



FIGURE 17 


OGIVE-CYLINDER SURFACE PRESSURE DISTRIBUTION 
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CIRCUMFERENTIAL ANGLE, 


TRAJECTORIES OF SEPARATION AND VORTEX CENTER POINTS 
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FIGURE 19 



COMPARISON OF NUMERICAL AND EXPERIMENTAL LIFT 
AND DRAG COEFFICIENT DISTRIBUTIONS 
OGIVE-CYLINDER BODY 
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FIGURE 20 
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VELOCITY VECTOR PLOT AT STA 1.135 



FIGURE 21 



VELOCITY VECTOR PLOT AT STA 10.0 



FIGURE 22 


VELOCITY VECTOR PLOT AT STA 12.0 
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FIGURE 23 



VELOCITY VECTOR PLOT AT STA 17.0 



FIGURE 24 



VELOCITY VECTOR PLOT AT STA 19.6 



FIGURE 25 



STATIC PRESSURE INSTRUMENTATION ON FUSELAGE FOREBODY 
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COMPARISON OF NUMERICAL AND EXPERIMENTAL 

SURFACE PRESSURES 

NUMERICAL RESULTS EXPERIMENT PLANE 



FUSELAGE STATION, in. 


FIGURE 27 



COMPARISON OF NUMERICAL RESULTS WITH EXPERIMENT 
CONTOURS OF LOCAL PITOT PRESSURE RATIO 
P Tp/ P Tco FUS STA 19.5 
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COMPARISON OF NUMERICAL RESULTS WITH EXPERIMENT 

CONTOURS OF LOCAL SIDESLIP ANGLE, /3 

FUS STA 19.50 




COMPARISON OF NUMERICAL RESULTS WITH EXPERIMENT 
CONTOURS OF LOCAL MACH NUMBER, M 

FUS STA 19.50 
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FIGURE 30 



COMPARISON OF NUMERICAL RESULTS WITH EXPERIMENT 
CONTOURS OF LOCAL ANGLE OF ATTACK, Q e 

FUS STA 19.50 




FIGURE 31 



CONTOURS OF LOCAL TOTAL PRESSURE RATIO, p s /p Sc0 

FUS STA 19.50 
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ASA-Lanelev. 1971 


LIFT AND DRAG COEFFICIENT DISTRIBUTIONS 

FUSELAGE CONFIGURATION 
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